Garnett构建自己的分类器以定义细胞类型

Original 周运来单细胞天地 2022-06-06

分享是一种态度

作者 | 周运来

男，

一个长大了才会遇到的帅哥，

稳健，潇洒，大方，靠谱。

一段生信缘，一棵技能树，

一枚大型测序工厂的螺丝钉，

一个随机森林中提灯觅食的津门旅客。

如果您的组织类型不存在分类器我们的仓库中，或者数据中不包含您期望的细胞类型，那么您需要生成自己的分类器。

训练分类器的第一步是加载单细胞数据。

因为Garnett 建立在 Monocle（http://cole-trapnell-lab.github.io/monocle-release）上，所以Garnett 的数据保存在CellDataSet (CDS)类的对象中。这个类派生自Bioconductor ExpressionSet类，它为那些分析过生物微阵列实验的人提供了一个常见的接口。Monocle提供了关于如何生成输入cds的详细文档here（http://cole-trapnell-lab.github.io/monocle-release/docs/#the-celldataset-clas）。

例如，Garnett包含一个来自PBMC 10x V1表达式数据的小数据集.

 1library(monocle3)
 2library(garnett)
 3# load in the data
 4# NOTE: the 'system.file' file name is only necessary to read in
 5# included package data
 6#
 7mat <- Matrix::readMM(system.file("extdata", "exprs_sparse.mtx", package = "garnett"))
 8fdata <- read.table(system.file("extdata", "fdata.txt", package = "garnett"))
 9pdata <- read.table(system.file("extdata", "pdata.txt", package = "garnett"),
10                    sep="\t")
11row.names(mat) <- row.names(fdata)
12colnames(mat) <- row.names(pdata)
13
14# create a new CDS object
15#pd <- new("AnnotatedDataFrame", data = pdata)
16#fd <- new("AnnotatedDataFrame", data = fdata)
17pbmc_cds <- new_cell_data_set(as(as.matrix(mat), 'sparseMatrix'),
18                              cell_metadata = pdata,
19                              gene_metadata = fdata)
20
21# generate size factors for normalization later
22#pbmc_cds <- estimateSizeFactors(pbmc_cds)#
23
24pbmc_cds = preprocess_cds(pbmc_cds, num_dim = 100)

构造标记(marker )文件

除了表达数据之外，您需要的第二个主要输入是一个标记文件(marker file)。标记文件包含以易于阅读的文本格式编写的单元类型定义列表。细胞类型定义告诉Garnett如何选择细胞来训练模型。每个细胞类型定义以“>”符号和细胞类型名称开头，后面是一系列带有定义信息的行。定义行以关键字和':'开头，条目用逗号分隔。

一个简单有效的例子:

1>B cells
2expressed: CD19, MS4A1
3
4>T cells
5expressed: CD3D

有几种方法可以在Garnett标记文件格式中定义细胞类型。通常，每个细胞的定义可以包含三个主要组件。只需要第一个组件。

细胞类型的第一个也是最重要的规范是它的表达式。Garnett提供了几种指定标记基因的选项，详情如下。

1Format  
2expressed: gene1, gene2 
3not expressed: gene1, gene2

这是为你的细胞类型指定标记基因的默认方法。在使用这个规范时，Garnett为每个细胞计算一个标记分数，包括潜在的泄漏（leakage）、总体表达水平和read 深度。expressed:标记应该特定于所定义的细胞类型。

Format	Example
expressed above: gene1 value, gene2 value	expressed above: MYOD1 4.2, MYH3 700
expressed below: gene1 value, gene2 value	expressed below: PAX6 20, PAX3 4
expressed between: gene1 value1 value2, gene2 value1 value2	expressed between: PAX6 10 20, PAX3 4 100

这是指定expression的另一种方法，如果您知道某个基因的确切表达范围，这种方法会很有用。但是，通常我们不建议使用这些规范，因为它们不会考虑每个细胞中的read深度和总体表达。数据值与输入数据的单位相同。

定义元数据

除了表达式信息之外，您还可以使用元数据进一步细化细胞类型定义。这也是指定数据中所期望的任何子类型的地方。

Format	Example
subtype of: celltype	subtype of: T cells
custom meta data: attribute1, attribute2	tissue: spleen, thymus

subtype of:允许您在定义文件中指定细胞类型是另一个细胞类型的子类型。
custom meta data:规范允许您为细胞类型提供任何进一步的元数据需求。CDS对象的pData表中的任何列都可以用作元数据规范。在上面的示例中，pData表中有一个名为“tissue”的列。

提供你的证据

最后，我们强烈建议您记录如何选择标记定义。为了便于跟踪，我们提供了一个额外的规范—references:—它将存储每种类型的引用信息。添加一组url或DOIs，它们将包含在分类器中。有关访问此信息的函数，请参见此处（https://cole-trapnell-lab.github.io/garnett/docs#viewing-references）。

添加注释

与R代码类似，我们已经包含了一个注释字符#，这样您就可以在您的标记文件中添加注释/注释。任何在同一条线上的#之后的内容都会被忽略。

所以一个完整的例子是这样的：

 1>B cells
 2expressed: CD19, MS4A1
 3expressed above: CD79A 10
 4references: https://www.abcam.com/primary-antibodies/b-cells-basic-immunophenotyping,
 510.3109/07420528.2013.775654
 6
 7>T cells
 8expressed: CD3D
 9sample: blood # A meta data specification
10
11>Helper T cells
12expressed: CD4
13subtype of: T cells
14references: https://www.ncbi.nlm.nih.gov/pubmed/?term=1200072

目前，标记文件不能有回车(\r)，如果您在某些Windows文本编辑器中生成标记文件，就会自动包括回车。相反，您必须使用换行字符(\n)。有许多方法可以在Windows上替换回车，例如使用notepad++，您可以使用replace、choose extended和replace all \r with \n。我打算尽快以更自动化的方式解决这个问题。

检查标记

因为定义标记文件通常是过程中最困难的部分，所以Garnett提供了一些函数来检查标记是否能够正常工作。相关的两个函数是check_marker和plot_marker。check_marker生成关于标记的信息表，plot_marker绘制最相关的信息。
除了包含的小数据集之外，我们还在包中包含了两个示例标记文件:pbmc_bad_markers.txt

在生成您的marker_check data.frame 之后，您可以使用我们内建的绘图函数plot_marker来查看结果。

1library(org.Hs.eg.db)
2marker_file_path <- system.file("extdata", "pbmc_bad_markers.txt",
3                                package = "garnett")
4marker_check <- check_markers(pbmc_cds, marker_file_path,
5                              db=org.Hs.eg.db,
6                              cds_gene_id_type = "SYMBOL",
7                              marker_file_gene_id_type ="SYMBOL")
8
9plot_markers(marker_check)

db: db is a required argument for a Bioconductor AnnotationDb-class package used for converting gene IDs. For example, for humans use org.Hs.eg.db（http://bioconductor.org/packages/3.8/data/annotation/html/org.Hs.eg.db）. See available packages at the Bioconductor website（http://bioconductor.org/packages/3.8/data/annotation/）. Load your chosen db using library(db). If your species does not have an AnnotationDb-class package, see here（https://cole-trapnell-lab.github.io/garnett/docs/#my-species-doesn-t-hav）.
cds_gene_id_type: This argument tells Garnett the format of the gene IDs in your CDS object. It should be one of the values in columns(db). The default is "ENSEMBL".
marker_file_gene_id_type: Similarly to above, this argument tells Garnett the format of the gene IDs in your marker file.

这个标记图提供了一些关于所选标记是否正确的关键信息。首先，红色标记“not in db”让我们知道标记ACTN在org.Hs.eg.db注释中没有作为“SYMBOL”出现。在本例中，它是一个打印错误。接下来，x轴显示每个标记的模糊度评分—当包含该标记时，测量有多少个cell接受了模糊标签—在本例中，ACTB和PTPRC具有很高的模糊度，应该排除。

check_marker输出的值和plot_marker绘制的值是分类器可以选择的cell 数量的估计值。然而，它使用启发式快速找到候选细胞，并不能完全匹配标记所选择的细胞。请使用这些数字作为相对的度量，而不是训练集的绝对表示。

关于歧义分数的进一步说明:歧义分数是当一个标记被包含在标记文件中时，一个标记被标记为“歧义”的cell 的分数。然而，一个高模糊度的分数并不一定意味着一个给定的标记是不具体的。这可能意味着一个不同的标记是罪魁祸首，但该标记也提名了许多其他未标记的细胞(高提名率)。在决定排除哪个标记之前，请仔细查看歧义度高的标记和最模糊的cell类型marker。

在制作标记图之后，您可能想要修改标记文件。在我们的示例中，我们将删除ACTN、ACTB和PTPRC以得到最终的“pbmc_test”。txt的标记文件。

默认情况下，Garnett 将基因id转换成ENSEMBL用于分类器。这使得分类器是可移植的，这样它们就可以对未来可能具有不同基因id的数据集进行分类。然而，对于某些生物体，ENSEMBL id要么不可用，要么不常用。在这些情况下，可以在train_cell_classifier和check_marker中使用参数classifier_gene_id_type来指定不同的ID类型。您选择的值将与分类器一起存储，因此在对未来的数据集进行分类时不需要再次指定它。

训练分类器

现在是训练分类器的时候了。参数应该与check_marker的参数非常接近。下面我将从默认值更改的一个参数是num_unknown参数。这告诉Garnett 应该比较多少个外群细胞。默认值是500，但是在这个只有很少cell的数据集中，我们需要更少的cell。

 1library(org.Hs.eg.db)
 2set.seed(260)
 3
 4marker_file_path <- system.file("extdata", "pbmc_test.txt",
 5                                package = "garnett")
 6pbmc_classifier <- train_cell_classifier(cds = pbmc_cds,
 7                                         marker_file = marker_file_path,
 8                                         db=org.Hs.eg.db,
 9                                         cds_gene_id_type = "SYMBOL",
10                                         num_unknown = 50,
11                                         marker_file_gene_id_type = "SYMBOL")
12head(pData(pbmc_cds))

 1> head(pData(pbmc_cds))
 2DataFrame with 6 rows and 4 columns
 3                           tsne_1           tsne_2
 4                        <numeric>        <numeric>
 5AAGCACTGCACACA-1  3.8403149909359 12.0841914129204
 6GGCTCACTGGTCTA-1 9.97096226657347 3.50539308651821
 7AGCACTGATATCTC-1 3.45952940410281 4.93527280576176
 8ACACGTGATATTCC-1 1.74394947394641 7.78267061846286
 9ATATGCCTCTGCAA-1 5.78344829514223 8.55889827553495
10TGACGAACCTATTC-1 10.7928530485958 10.5852739146963
11                       Size_Factor FACS_type
12                         <numeric>  <factor>
13AAGCACTGCACACA-1 0.559181445161514   B cells
14GGCTCACTGGTCTA-1 0.515934033527584   B cells
15AGCACTGATATCTC-1 0.698028398302026   B cells
16ACACGTGATATTCC-1 0.815631008885519   B cells
17ATATGCCTCTGCAA-1  1.11532798424345   B cells
18TGACGAACCTATTC-1 0.649469901028841   B cells

运行train_cell_classifier之后，类型为“garnett_classifier”的输出对象包含对细胞进行分类所需的所有信息。

查看分类基因

Garnett 分类是使用多项弹性-网络回归训练（multinomial elastic-net regression）。这意味着选择某些基因作为区分细胞类型的相关基因。所选择的基因可能是有趣的，所以Garnett 包含了一个访问所选择基因的功能。注意:Garnett 没有对输入标记进行正则化，所以无论如何，它们都会被包含在分类器中。

我们用来查看相关基因的函数是get_feature_genes。参数是分类器，您想查看哪个节点(如果您的树是分层的)—使用“root”作为顶部节点，使用父细胞类型名称作为其他节点，使用db作为您的物种。如果设置convert_ids = TRUE，该函数将自动将基因id转换为符号。

 1feature_genes <- get_feature_genes(pbmc_classifier,
 2                                   node = "root",
 3                                   db = org.Hs.eg.db)
 4head(feature_genes)
 5                   B cells      T cells     Unknown
 6(Intercept)     -4.44471924  2.590206989  1.85451225
 7ENSG00000187608 -0.10903325 -0.199658082  0.30869133
 8ENSG00000157873  0.06411921 -0.295797377  0.23167817
 9ENSG00000157191 -0.07677221 -0.155893113  0.23266532
10ENSG00000173436  0.02649831  0.089393639 -0.11589195
11ENSG00000117318  0.07400983 -0.009415156 -0.06459467
12

Viewing references

我们在上面解释了如何在标记文件中包含关于如何选择标记的文档。为了获取这些信息—查看如何为已经训练好的分类器选择标记—使用函数get_classifier_references。除了分类器之外，还有一个额外的可选参数，称为cell_type。如果传递细胞类型的名称，则仅打印该单元类型的引用，否则将全部打印。

1get_classifier_references(pbmc_classifier )
2$`B cells`
3[1] "https://www.ncbi.nlm.nih.gov/pubmed/?term=23688120"
4[2] "https://www.ncbi.nlm.nih.gov/pubmed/?term=21149806"
5
6$`T cells`
7[1] "https://www.ncbi.nlm.nih.gov/pubmed/?term=1534551"